home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_2 / steff.m < prev   
Encoding:
Text File  |  1994-06-05  |  1.5 KB  |  63 lines  |  [MATF/MATL]
  1. function [p,yp,err,Q] = steff(f,df,p0,delta,epsilon,max1)
  2. % [p,yp,err] = steff(f,df,p0,delta,epsilon,max1)
  3. % [p,yp,err,Q] = steff(f,df,p0,delta,epsilon,max1)
  4. % Steffensen's method is used to locate a root.
  5. % f is the function, input.
  6. % d is the derivative of f, input.
  7. % p0 is the starting point, input.
  8. % delta is the tolerance for p, input.
  9. % epsilon is the tolerance for yp, input.
  10. % max1 is the maximum number of iterations, input.
  11. % p  is the root, output.
  12. % yp is the function value, output.
  13. % err is the error estimate for p, output.
  14. % Q is the is the vector of iterations, output.
  15. y0 = feval(f,p0);
  16. p = p0; yp = y0;
  17. Q(1) = p0;
  18. for k=1:max1,
  19.   df0 = feval(df,p0);
  20.   if df0 == 0,
  21.     dp = 0;
  22.   else
  23.     dp = y0/df0;
  24.   end
  25.   p1 = p0 - dp;
  26.   y1 = feval(f,p1);
  27.   p = p1; yp = y1;
  28.   Q = [Q,p1];
  29.   err = abs(dp);
  30.   relerr = err/(abs(p1)+eps);
  31.   if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
  32.   df1 = feval(df,p1);
  33.   if df1 == 0,
  34.     dp = 0;
  35.   else
  36.     dp = y1/df1;
  37.   end
  38.   p2 = p1 - dp;
  39.   y2 = feval(f,p2);
  40.   p = p2; yp = y2;
  41.   Q = [Q,p2];
  42.   err = abs(dp);
  43.   relerr = err/(abs(p1)+eps);
  44.   if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
  45.   d1 = (p1 - p0)^2;
  46.   d2 = p2 - 2*p1 + p0;
  47.   if  d2 ==0, 
  48.     dp = p2 - p1;
  49.     p3 = p2;
  50.     break;
  51.   else
  52.     dp = d1/d2;
  53.     p3 = p0 - dp;
  54.   end
  55.   p0 = p3;
  56.   y0 = feval(f,p0);
  57.   p = p0; yp = y0;
  58.   Q = [Q,p0];
  59.   err = abs(dp);
  60.   relerr = err/(abs(p3)+eps);
  61.   if (err<delta)|(relerr<delta)|(abs(y1)<epsilon), break, end
  62. end
  63.